Vortices, zero modes and fractionalization in bilayer-graphene exciton condensate 
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A real-space formulation is given for the recently discussed exciton condensate in a symmetrically 
biased graphene bilayer. We show that in the continuum limit an oddly-quantized vortex in this 
condensate binds exactly one zero mode per valley index of the bilayer. In the full lattice model 
the zero modes are split slightly due to intervalley mixing. We support these results by an exact 
numerical diagonalization of the lattice Hamiltonian. We also discuss the effect of the zero modes 
on the charge content of these vortices and deduce some of their interesting properties. 



Introduction. — Nontrivial topological configurations 
play a significant role in our understanding of a wide 
range of physical systems and often have interesting and 
potentially useful properties. In one spatial dimension, 
solitons and instantons interacting with fermions often 
induce fractional fermion numbers [TJ [2] as manifested 
in the canonical example of the polyacetylene chain [3] . 
In two dimensions, an important example is provided by 
vortices in a superfluid or a superconductor, precisely 
quantized to make possible the coexistence of the conden- 
sate with supercritical macroscopic rotation or magnetic 
field, respectively. Such vortices may carry fermonic zero 
modes which typically imply a host of unusual properties. 
In a planar spin-polarized p + ip superconductor half- 
quantum vortices can exist and have been argued to obey 
non-Abelian exchange statistics due to the presence of an 
unpaired Majorana fermion bound to the vortex 0] [S]. 
More recently similar physics has been shown to occur in 
certain simple models describing fermions on the honey- 
comb and TT-flux square lattices where vortices in a bond 
dimerization pattern carry fractional charge |BJ |7] and 
obey Abelian fractional statistics O |9] . 

Topological structures in the above 2D examples and 
their associated phenomena afford a special degree of pro- 
tection against local perturbations and have been iden- 
tified as the key resource for fault-tolerant quantum in- 
formation processing |10j . It would be clearly desirable 
if such systems could be produced and studied in the 
laboratory. Unfortunately, prospects for the experimen- 
tal realization of the above systems are not particularly 
good: no spin polarized p+ip superconductors are known 
to exist and it is not quite clear how one would realize 
the dimerization patterns required in Refs. [6,7. 

In this Letter we discuss a physical system that is likely 
to be produced in a laboratory in the near future and 
show that it exhibits some of the key features of the 
aforementioned 2D models [3 |5J [71 [HI [S] . Specifically, we 
consider a graphene bilayer separated by a dielectric bar- 
rier yjj depicted in Fig. [ija). When biased by external 
gates the perfect nesting of the electron and hole Fermi 
surfaces in different layers creates ideal conditions for the 
formation of an exciton condensate (EC) i.e. the coherent 
hquid of bound states of electrons in one layer and holes 
in the other. It has been argued recently that for realis- 
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FIG. 1: (Color online) (a) The schematic structure of the 
system, and the proposal to create and effective axial mag- 
netic field, (b) The spectrum near the nodal point where the 
dashed (green) line shows the non-interacting, and the solid 
(red) line the interacting case with the shifted nodes in exter- 
nal potential (V) and the exciton gap (m) marked. 



tic values of parameters the EC can exist even at room 
temperature and that it depends only weakly on the 
exact stacking of the two layers [12]. In what follows we 
construct a real-space model for this system and use it 
to study the internal structure of a vortex in the com- 
plex order parameter characterizing the EC. We show 
that it contains two fermion zero modes, one for each 
valley index of the bilayer, which are split slightly due 
to inter-valley mixing. Although there is no net charge 
associated with a vortex at half filling, the vortex binds 
fractional "axial charge" , defined as the difference of the 
charge between the layers. We argue that such a vortex 
will obey fractional exchange statistics and comment on 
possible experimental signatures of our findings. 

Exciton condensate. — We consider the simplest model 
containing the essential physics of the system at zero tem- 
perature. It consists of fermions hopping between nearest 
neighbor (n.n.) sites on each of two layers of honeycomb 
lattice stacked directly on top of each other. Since the 
electron spin plays no role in the formation of EC we 
consider only a single spin projection. There is no direct 
hopping between the two layers. EC formation is driven 
by the interlayer Coulomb repulsion which we model by 
an effective short-range repulsion between the fermions 
in different layers at the same planar site. As argued pre- 
voiusly [TP a layer of graphene is in the semi-metallic 
phase, meaning that the intra-layer Coulomb interaction 
is below the critical value needed for charge ordering. 
The long-range part of the interlayer Coulomb repulsion 
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should only enhance the formation of the EC. 

The Hamiltonian is, thus, H = Hi + H2 + Hu, 
where the in-plane Hamiltonian Ha — —tj^/ 



{—)°'VJ2i''^ia for layers a = 1,2 and the interaction 
Hu — U '^^niini2. Here Cia annihilates an electron at 

site i and layer a and riia = cl^c^^- 

The mean-field order parameter of the EC is = 
U{cl2C^l). Decoupling Hjj in this channel we find 
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HmF = Hi + H2-J2 (^»4lCi2 + h 
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We now define A± ~ ^{^A ± As) where A and B are 
the two sublattice sites in the same unit cell, and the 
spinor field operator ipi — {cBii,—CAii,CBi2,CAi2)'^- For 
a uniform A± we may write the Hamiltonian compactly 
in momentum space as H = V'k^kV'k + ^Oi with Eq 
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|A_| ), N the number of sites per layer. 



hk = 70 



7iRetk + 72lnitk + V^7o75 + 



-i7i72|A4.|e^ 



175 x+ 



(2) 



Here ^k = -^Eaen.n. 
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are Dirac matrices in the Weyl representation, 75 = 
— i7o7i7273 = 1T3 §5 1, and we have used a polar rep- 
resentation A± = |A±|e'^±. 

At nonzero A_ a finite energy gap. Eg < |A_|, opens 
at half filling with the maximum occurring for fixed A_ 
when |A^| — \Ab\- The lowest ground state energy oc- 
curs when the gap is maximum. Since a finite |A+| does 
not open a gap on its own we also expect the ordered 
state with A_|_ = to set in first. This can be checked 
explicitly by a calculation of the critical interaction for 
|A+|, which, for small V < yields Uc+ /V > t. 

For A_|_ = 0, the spectrum is ±i?ka, 



i?k. = V(l^k|-(-W + l^-P- (3) 
This is shown in Fig. [ijj. The gap equation reads 27V = 
U X^ka ^ka where the sum is over the occupied states. 
The critical value of the interaction is given by 2N /Uc- — 
EL IKkl - Around the nodes the right-hand 

side diverges logarithmically for a — 2. So, Uc- = as 
expected, in contrast to Uc+- The dependence of |A_| 
on U and V can be found in the nodal approximation to 
be 
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with an ultraviolet cut-off A ~ t and in the limit | A_ | <C 
V <^ k <^t^ /U. The state with A+ = and A_ 7^ is 
the EC phase considered in Refs. [TTJ[T2]. 

Low-energy theory. — At long wave-lengths the system 
is dominated by the excitations around the two inequiva- 
lent valleys K± = (±47r/3\/3a, 0) of the graphene layers. 




where tK±+p ~ '^Px+Wy (We have set the Fermi veloc- 
ity 3ta/2 =1.) The linearized mean-field Hamiltonian for 
these excitations is h^ + h_ with = 7173/1-7371 = Ti, 

W = 70 {jiPx + l2Py + 1^7075 + |m|e"'''^^) , (5) 

where p = — iV is the momentum operator and we have 
replaced A_ = m— Imje'-*^ for ease of use. The mass m is 
now also taken to vary with position. Since the valleys are 
decoupled in the low-energy theory we shall first study a 
single valley. We discuss the effect of intervalley mixing 
later and return to the full Hamiltonian in the numerics. 

The spectrum of H is symmetric around zero. This is 
seen by noting that 

72^*72 = n, (6) 

and 7I = —1. Thus, for every eigenstate ipE oi energy 
E the state ^ipE = 72^^; is an eigenstate of energy —E. 
This symmetry is furnished by the antiunitary operator 
J7 = fjt = ^^if where K is the complex-conjugation op- 
erator: {ri,7i} = 0. When V = Q formally, there is 
another unitary operator, 7073, that anticommutes with 
Ti (Ref. [Ojl . For V 7^ 0, in contrast, Q. is the only an- 
ticummuting operator. 

Zero modes. — We will now consider a vortex in the 
EC order parameter, i.e. m{r — s- 00, 9) = woe'"^ where 
n G Z is the vorticity, and show that, within the low- 
energy theory, there is exactly one zero mode per valley 
bound to the vortex for odd n. 

Before presenting the explicit solution, we note the fol- 
lowing. For V = Eq. ^ formally coincides with the 
Hamiltonian studied in Refs. j6l[7l[T3], and hence there 
are |n| exact zero modes for n-fold vortex [13]. We now 
imagine slowly turning on V and note that far from the 
vortex center the spectrum must remain gapped. Ow- 
ing to Eq. (|6| , which requires a symmetric spectrum, we 
could at most split an even number of zero modes. There- 
fore, at least one zero mode must survive for odd n. Also, 
since the spectral symmetry for V ^ is generated by an 
antiunitary operator, similar to the situation with half- 
quantum vortices in a, p + ip superconductor, we do not 
expect more than one zero mode to survive [TS]. This is 
confirmed by our explicit solution below. 

To find zero-mode solutions Hipo = we limit our- 
selves for simplicity to the quantum limit, m(r ^ 0,0) — 
mge'"^. The case of a more general radially symmetric 
vortex could be treated similarly. In the zero-energy sub- 
space, [H, il] = 0. As a result we may take ^ipo = ^''Po 
where the phase A = 1/A* is an "eigenvalue" of antiuni- 
tary f2. We can always choose such a^'o' if ^V'o 7^ -^V'o we 
can pick tp^ = [I + \*n)^o for which QtpQ = XtpQ. Since 
for an "eigenstate" ^(t)\ — \4>\ we have ri(e~'^'^0A) — 
(e''^A)(e"'''/^0A) we may further limit ourselves to A = 1. 
Such cigcnstates are of the form ■00 = (/)5,iff*7 — i/*)^ 
and form an over-complete basis of the spinor space. 

With this preparation we find two independent equa- 
tions for our zero modes, 
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Any constant phase e'^ of mo can be absorbed by 
if, 9) e'l^' {f,9) so we take < mo € R. We 
now make a "one-phase" ansatz f{r,6) = F{r)e^°-^ and 
g{r,9) = G{r)e^^^. We have checked that a "two-phase" 
ansatz like the one used in Ref. 14 gives no new solutions 
when V ^ 0. Thus, we are left with at most a single 
zero mode. We may eliminate the angular dependence 
by taking a = b — 1 — . We see that when n is even 
there is no single-valued solution. The radial part yields 
the bound-state Bessel-function solution 
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For V and n = ±1 we recover the zero mode of 
Ref. [6] which has support only on one sublattice in each 
layer. Since we have at most a single zero mode, the 

V ^ limit when |n| > 1 is singular. 

A localized zero-mode in a gapped, particle-hole sym- 
metric system is known to carry a fractional charge ±e/2 
[TJ 121 |3J [SJ [J . Since we have two independent valleys, 
however, we expect two zero modes to exist for odd n. 
Thus, charge fractionalizion should not occur. Neverthe- 
less, we show below that the zero modes, which are in 
addition split slightly due to intervalley mixing, lead to 
fractional "axial charge" bound to a vortex. This axial 
charge is experimentally measurable and also endows the 
vortex with fractional exchange statistics. In the real sys- 
tem we must also include the spin of electrons. However, 
unlike the valley degree of freedom the vortex need not 
mix spins since in principle we can have a vortex in just 
one spin species. 

Numerics. — We have performed exact diagonalization 
of the mean-field Hamiltonian ([ij on honeycomb lattices 
with up to 51 X 30 sites per layer. (A layer with x 

Ny sites has spatial dimensions ^N^a x '^NyU.) Our 
numerical results support the conclusion that near-zero 
modes exist only for odd vorticity. 

At half filling the charge is balanced between the lay- 
ers, so that the total charge is always uniform, with or 
without a vortex. However, the charge difference, 6Q, be- 
tween the layers, which we call the "axial charge" in anal- 
ogy with Ref. |9l shows interesting features. (The axial 
charge is proportional to the electric dipole between the 
layers.) In the uniform system and at mo = we have the 
axial density Sp = SQ/N = j^E\t,\<v'^ « ^^(^A)' 
for y <C Fig. [2] shows our numerical results for the 
axial charge of a typical system with 39 x 22 x 2 sites. 
In the uniform system, we observe a nearly linear depen- 
dence of (5/5 on V < mo which turns into quadratic for 

V > niQ. The density profile in Fig. ^c) shows oscilla- 
tions around the edges and dependence on the type of 
the edge (zig-zag for the horizontal and arm-chair for the 
vertical edges). 

A vortex sitting at the center of the system binds an 
irrational fraction of the axial charge, SQy, on top of the 
uniform background. The dependence of SQ^ on V and 
mo in Fig.|2jb) shows that the value as well as the sign of 
this extra axial charge may be selected continuously by 
tuning external parameters. The density profile of this 
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FIG. 2: (Color online) The "axial charge" . (a) Density in 
the uniform system, and (b) the extra charge bound to a 
single vortex as a function of V and mo. The corresponding 
density profile for V/t = 0.4 and mo/t = 0.3 in (c) the uniform 
system, and (d) the system with a vortex at the center. In 
all plots, each layer has 39 x 22 sites with a lattice size ~ 
(35)^. The lattice sites in (c) and (d) are at the center of the 
triangles. 



extra axial charge follows the shape of the zero-mode 
wavefunction and is displayed in Fig. [2|d). 

Fig. |3] shows the energy splitting of the zero modes 
obtained numerically. We can understand this analyti- 
cally by evaluating the matrix element between the zero- 
modes obtained in the single-valley approximation and 
then diagonalizing the resulting 2x2 Hamiltonian in the 
zero-mode subspace. This yields the zero mode splitting 
zblel with 
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/(l/P + l5P)dr ' 



(9) 



where Q is the inter- valley momentum. Eq. (9| shows the 
same qualitative behavior as the numerical results. In 
order to obtain a reasonable quantitative fit we must re- 
place Q with an "effective" internodal wavevector, jy^^Q, 
which represents the effects of the full spectrum. (Note 
that Q is completely extraneous to the nodal approxi- 
mation.) A good fit is obtained with 77 w 3.0 for low to 
intermediate values of mo and V. 

Fractional statistics. — The vortex carries fractional ax- 
ial charge 6Qy and, upon performing a singular gauge 
transformation employed in Ref. 9, it can be seen to carry 
a half of an axial flux quantum, defined as the differ- 
ence between the effective magnetic fluxes piercing the 
two layers. Employing the standard argument |17| such 
charge-flux composite particles are expected to behave as 
Abelian anyons with the exchange phase 7 = ttSQ^. In 
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FIG. 3: (Color online) (a) The zero-mode energy splitting. 
The squares show our numerical results for a system size ~ 
(35)^. The sohd lines show the fitting with Eq. (b) 
The exchange phase vs. axial charge. The solid (red) line 
is 7 = ttSQv 



order to verify this hypothesis we have computed the ex- 
change phase numerically using a generalized Bargmann 
invariant jl6l. Our results, summarized in Fig. [3'b), are 
consistent with the expected exchange phase, although 
system size limitations result in relatively large error 
bars. It is important to emphasize that unlike the sit- 
uation in the fractional quantum Hall fluids the value of 
the fractional charge and the exchange phase here is not 
protected by symmetry or topology but depends contin- 
uously on externally adjustable parameters. In this sense 
our situation is similar to "irrational fractionalization" of 
charge and statistics discussed in Ref. 18] 

Discussion. — The exciton condensate considered in 
this work should exhibit a host of unusual behaviors ob- 
servable experimentally. The system is an insulator in 
the total charge channel but a superconductor in the ax- 
ial channel. This behavior should lend itself to direct ob- 
servation in transport if the two layers could be indepen- 
dently contacted. Such measurements are possible in bi- 
layers fabricated using semiconductor structures [THl [TH] . 



The exciton order parameter has been shown to be 
minimally coupled to the "axial gauge field" i.e. the 
part of the vector potential A that produces the anti- 
symmetric combination SB =(Bi — B2)_Lof the mag- 
netic field normal to the plane of the layers 130]. Vor- 
tices in the EC therefore may be generated by the axial 
field just like Abrikosov vortices are generated by ordi- 
nary magnetic field in a superconductor. Indeed it is 
possible to arrange the external magnetic field in a way 
that 6B is much greater than the symmetric combina- 
tion (Bi + B2)_L. One such arrangement is schematically 
shown in Fig. [ija). Using standard arguments unidirec- 
tional drift of vortices through the system can be seen to 
produce an axial Hall voltage from which the fractional 
axial charge carried by an individual vortex can be in- 
ferred. Such a drift can occur in response to the Magnus 
force acting on vortices in external applied axial current 
or, alternatively, due to a temperature gradient applied 
across the sample. 

The zero modes found here exist in the excitonic in- 
stability of a Fermi surface. The existence of the under- 
lying Fermi surface is at the root of excitonic instability 
at infinitesimal coupling and also causes the oscillatory 
Bessel-function behavior of our solution. Previous zero- 
mode solutions have, to the best of our knowledge, only 
appeared in systems with Dirac points. It is very in- 
teresting to see if zero modes of this type can be found 
in other physical systems. We also note that the U(l) 
vortices considered here are allowed by lattice symme- 
tries, hence the issue of linear confinement occuring in 
Refs. [6, 7 does not arise. Yet, vortices in the present 
model exhibit some of the same fascinating properties. 
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